Density Profiles in Open Superdiffusive Systems 
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We numerically solve a discretized model of Levy random walks on a finite one-dimensional 
domain in the presence of sources and with a reflection coefficient r. At the domain boundaries, the 
steady-state density profile is non-analytic. The meniscus exponent fi, introduced to characterize 
this singular behavior, uniquely identifies the whole profile. Numerical data suggest that /i = 
a/2 + r(a/2 — 1), where a is the Levy exponent of the step-length distribution. As an application, 
we show that this model reproduces the temperature profiles obtained for a chain of oscillators 
displaying anomalous heat conduction. Remarkably, the case of free-boundary conditions in the 
chain correspond to a Levy walk with negative reflection coefficient. 
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Anomalous diffusion and transport are ubiquitous phe- 
nomena that occur in many complex systems, ranging 
from molecular motion to human mobility In the 
framework of random walk theory much progress has 
been made in the last years in the characterization of 
superdiffusive motion in infinite domains [l], 0| • Exten- 
sion to the case of finite systems is relatively less devel- 
oped 0, Typically, only the simple case of absorb- 
ing Boundary Conditions (BC) without sources has been 
treated @ . An even less studied problem is that of fi- 
nite and open superdiffusive systems where particles and 
energy arc steadily exchanged with some external reser- 
voirs [6|]. This setup describes a variety of systems with 
anomalous transport in the off-equilibrium regime. For 
many-particle systems this is exemplified by anomalous 
heat conduction in low-dimensional lattices 0, @| , as well 
as by enhanced thermodiffusion effects 0. Moreover, 
there exist at least two physical setups where the the- 
oretical predictions can be experimentally tested: light 
superdiffusion 1^, 11 1 and anomalous heat conduction in 



individual nanotubes 12 1 



In the present Letter we study the spatial distribution 
of the relevant field (particle density, temperature, light 
intensity, depending on the context) in the presence of an 
external gradient. By numerically investigating a model 
of particles performing Levy flights or walks, we find that 
the shape of the profile is nonlinear in the thermody- 
namic/continuum limit. This is due to the long-range 
correlations which affect the transport far away from the 
boundaries. A comparison with the temperature profiles 
obtained in a model of heat conductivity indicates that 
the shape is fully determined by the exponent of the un- 
derlying Levy process and by the BC. However, the con- 
nection between the (free/fixed) BC of the conduction 
problem and those of the diffusive process (reflectivity) 
is not trivial. We indeed find that free BC correspond to 
a negative reflectivity. 

The model - Let i denote the position of a discrete-time 
random walker on a finite one-dimensional lattice (1 < 



i < N). In between consecutive scattering events, the 
particle either jumps instantaneously (Levy flight - LF) 
or moves with unit velocity (Levy walk - LW |l3|) over a 
distance of m sites, that is randomly selected according 
to the step-length distribution 
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^0 = 0. 
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Here a is the Levy exponent and q is a normalization 
constant. We limit ourselves to the case 1 < a < 2, where 
the averages are finite but the variance not. We formulate 
the problem by introducing the vector W = {Wj(t)}, 
where W% is the probability for the walker to undergo a 
scattering event at site i and time t. It satisfies a master 
equation, which, for LFs, writes 



W(i + l) = QW(i) + S, 



(2) 



where S accounts for the particles steadily injected from 
external resorvoirs (see Fig. [T] for a pictorial representa- 
tion) ; Q is a matrix describing the probability of paths 
connecting pair of sites. In the simple case of absorbing 
BC, it is readily seen that Qji is equal to the probability 
ipj-i of a direct flight, as from Eq. ^ (see Fig.[T]). In the 
LW case, the W components in the r.h.s. must be esti- 
mated at different times (depending on the length of the 
path followed from j to i). Since, the stationary solution 
is the same in both setups, this difference is immaterial, 
and we thereby refer to LFs, since Eq. ([2]) is easier to 
solve iteratively. 

In the LF case, Wj is equal to the density P$ of particle 
at site i. This is no longer true for the LW, as includes 
those particles that are crossing during a ballistic step. 
In the simple case of absorbing BC, one can write 



Pi = Wi 



N 

i=i 



(3) 



account for the 



where Ri = J2j>i V'i an d h = Y^j 

particles which, having started their walk respectively 
inside or outside the domain, transit at the zth site. 
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FIG. 1: (Color online) Scheme of the model: an ensemble 
of Levy walkers on a finite one-dimensional lattice interact- 
ing with one external reservoir that inject particles with a 
prescribed constant rate Si at each site. 



The source term is fixed by assuming that the reservoir 
is a semi-infinite lattice, homogeneously filled by Levy 
walkers of the same type as those residing in the domain. 
This amounts to defining S m — sm~ a , where s mea- 
sures the density of particles and m the distance from 
the reservoir. It is easy to verify that in the presence of 
two identical reservoirs at the lattice ends, the density is 
constant (for any TV), showing that our definition satisfies 
a kind of "zeroth principle" as it should. 

In the nonequilibrium case, it is not necessary to deal 
with two reservoirs. The linearity of the problem teaches 
us that it is sufficient to study the case of a single reser- 
voir, that we assume to be in i — 0: the effect of, say, a 
second one on the opposite side can be accounted for by 
a suitable linear transformation. For values of N large 
enough, the steady-state density depends on i and N 
through the combined variable x = i/N, Pi = P(x). As 
seen in Fig. [2 P vanishes for x — »• 1 because on that side 
the absorbing boundary is not accompanied by an incom- 
ing flux of particles. Besides, we have also verified that 
Wi oc Pi in agreement with the expectation that LFs and 
LWs should have by the same stationary properties. 

Fractional Laplacian- For a more accurate consistency 
check, we now compare the solution of our model with 
that of the stationary Fractional Diffusion Equation 
(FDE) D2P = -ct(x) on the interval < x < 1 (see 
e.g. Ref. [14| and references therein for a precise defini- 
tion of the operator D%). The source term a(x) must 
be chosen so as to describe the effect of the external 
reservoirs. A condition to be fulfilled is that two identi- 
cal reservoirs yield a homogeneous state P(x) = const. . 
A straightforward calculation using the integral defini- 
tion of the operator 14J shows that this happens for 



a(x) 



,(*) 



+ (1 — x) a (we, henceforth, ig- 



nore irrelevant proportionality constants) 22]. It is thus 
natural to associate o~(x) — x~ a to the nonequilibrium 
case with a single source in x — 0. The numerical solution 
of the FDE is reported in Fig. [2] (dashed line). It over- 
laps almost perfectly with the stationary solution of the 
discrete model, thereby confirming the correctness of our 
choice. Moreover, this shows that long-ranged sources 
are needed to reproduce the superdiffusive profiles in the 
continuum limit. 
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FIG. 2: (Color online) Density profiles as a function of the 
variable x = i/N for absorbing BC and a = 3/2. Solid line: 
solution of the Eq. © with N = 4000. Dashed line: solution 
of the FDE discretized on a mesh of M = 800 points; the 
latter curve is rescaled by a suitable constant. The inset is a 
magnification of the leftmost region. 



Boundary singularity - A distinctive feature of the pro- 
file is that it is not analytic at the boundaries. Indeed, 
the data for x — > are well fitted by 



P(x) = P(0) + Cx>* 



(4) 



(the same behavior occurs for x — > 1, as the profiles are 
symmetric). In view of the similarity with the shape of 
the liquid surface close to a wall, we metaphorically term 
fi as the meniscus exponent. The density drop is akin to 
the well-known Kapitza resistance, which has been pre- 
viously discussed for one-dimensional systems [l5[ . How- 
ever, the nonanalytic behavior is peculiar of anomalous 
kinetics, as opposed to the familiar linear shape in stan- 
dard diffusion. For the above discussed case of absorbing 
boundaries, we find that [i ~ a/2. This value is consis- 
tent with the singular behavior of the eigenfunctions of 
the fractional Laplacian 14 1 . 

Reflecting boundaries - We now generalize the discrete 
model, by assuming that any time a particle reaches the 
boundary, it is reflected with probability r (0 < r < 1). 
This extension is of central interest for real problems, e.g. 
in anomalous diffusion of light but we have to pay 
the price of a more complicate mathematical treatment 
to account for multiple reflections. In fact, the transi- 
tion probability Qji from the site j to site i is the sum 
of infinitely many terms that correspond to paths un- 
dergoing multiple reflections. By calling B>i and B r the 
left and right boundaries, respectively, and viewing any 
path as a sequence of symbols that identify the start- 
ing and final site, as well as all reflection points, the 
four families that contribute to Qji can be encoded as 
(for 1 < j < i < N - 1): j(B r Bi) m i, jBi{B r Bi) m i, 
jB r (BiB r ) m i, and j(BiB r ) m+1 i with m = 0, 1, . . .. Sym- 
metric relations hold for j > i, while special expressions 
must be invoked for the degenerate cases (i,j = 1, N — 1, 
i = j). The weight of a path of length I is then given by 
r k ipi, where k is the number of reflections (k = 2m + 1 
in the second of the four families). By summing up all 
contributions for the same pairs of sites, we obtain Qji. 
A similar treatment must be applied to the source term 
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to take into account that each particle entering from ei- 
ther side may be trapped for an arbitrarily long time. 
Another complication is that for LWs, the expression for 
P{x) requires including all paths arising from multiple 
reflections. Notice however, that all such intricacies do 
not change the rank of the matrix Q. 

Unexpectedly, we found that the shape of the profile 
depends on r for fixed a. This is entirely new with re- 
spect to normal diffusion, when the profile is linear irre- 
spective of the boundary conditions. The meniscus expo- 
nent is a convenient tool to parametrize this functional 
dependence on r. The estimated numerical values of \x 
for a — 3/2 are reported in Fig. [3j Once more, we found 
that Wi oc Pi and that, accordingly, both profiles give 
compatible values of \i (see symbols in Fig. [3]). By as- 
suming a linear dependence of \i on both r, and a, we 
conjecture 

'a 
.2 



advantage of allowing for an exact solution of the associ- 



ated Fokker-Planck equation [21 



i« = 



1 



(5) 



This expression is consistent with the \x = a/2 value 
found above for r = 0. Moreover, for a = 2 (normal 
diffusion) it yields ji = 1, as it should. The discrepancy 
between formula (|5|) and the numerical values is every- 
where compatible with the finite-size corrections (data 
has been obtained with N = 8000). 
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FIG. 3: (Color online) Dependence of the meniscus exponent 
fi on the reflection probability r for a = 3/2. Full circles and 
stars are measures from fitting of P{x) and W(x), respectively 
(see text). The dashed line is given by the formula ©. 

Application to the heat conduction problem- Tem- 
perature profiles in one-dimensional systems displaying 
anomalous energy transport [§] represent a good test- 
ing ground of the above approach. In fact, energy per- 
turbations do spread like in a Levy process [17j . Also, 
for an infinite chain of harmonic oscillators with energy 
and momentum conserving noise [18| , it was proven that 
the phonon Boltzmann equation reduces to a FDE with 
a = 3/2 [l9[ . We thus expect model J2) to apply to a 
finite chain of oscillators, (i) coupled with two Langevin 
heat baths (with a damping constant A), and (ii) with 
random collisions that exchange the velocities of neigh- 
boring particles with a rate t [20|| . This model has the 




FIG. 4: (Color online) Temperature profile of the oscillator 
chain with conservative noise with free boundary condition 
and A = 7 = 1 (solid line) and solution of the master equation 
with reflection coefficient r = —0.1 (dashed line). 

In Fig. U we compare the temperature profile T(x) 
(suitably shifted and rescaled) of the heat-conduction 
model with free BC and the solution of our discrete Levy 
model with a reflection coefficient r = —0.1. Since they 
are essentially indistinguishable, we can conclude that 
the Levy interpretation does not only allow to explain 
the anomalous scaling of heat conductivity [ljj , but also 
the peculiar shape of T{x). The weird (negative) value 
of r requires a specific discussion. Since Q is the sum of 
terms proportional to r k (k being the number of reflec- 
tions from the boundaries), even and odd A:s yield, for 
r < 0, positive and negative contributions, respectively. 
Thus, we let Q = Q e Q . The presence of an unphysical 
negative term can be avoided by introducing two families 
of walkers (W + and W _ ) and interpreting the reflection 
as a change of family. Under these assumptions, the evo- 
lution equations read 



W±(t + 1) = Q e W ± (i) + Q„W T (t) + S d 



(6) 



where all terms are positive as they should. Under the 
further assumption that S + = S and S~ = 0, we can 
easily see that W = W + — W satisfies the original 
Eq. |(3J). In other words, the relevant quantity to look at 
is the difference between the densities of the two different 
families. Why it is necessary to invoke the presence of 
such two families and what is their physical meaning in 
the context of heat conductivity is an open problem. 

In the case of a chain with fixed BC, the temperature 
profile T(x) can be computed analytically 2l| and it is 
thereby found that fi = 1/2. By inserting this value in 
Eq. ([5|) and recalling that a — 3/2, we find that r = 1, 
i.e. the fixed-BC T{x) corresponds to the case of per- 
fectly reflecting barriers. Unfortunately, this (physically 
reasonable) result could not be tested quantitatively. In- 
deed, it turns out that finite-size corrections become in- 
creasingly important upon increasing r, and for r close 
to 1, it is practically impossible to achieve convergence 
to the steady-state. 

Local equilibrium- In the context of LWs, we can go be- 
yond the determination of the shape of P and check the 
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existence of local equilibrium. To this aim, we decompose 
Pi as J2 n [Pi R ( n ) + Pt{ n )\ where P^' R {n) is the proba- 
bility that a walker on site i and moving respectively to 
the left or right will undergo the next scattering event 
after n time steps. Then, we compare P^' R (n) with its 
equilibrium value P eq (n), obtained from the stationary 
solution of Eq. ([2]) with two equal sources on both sides 
[P eq (n) is independent of i and of L, R}. Unsurprisingly, 
P^' R {n) ~ P eq {n) close to the left boundary. Moving 
away from the reservoir, the densities Pi{n) start to de- 
viate from P eq (n) (see respectively the broken and solid 
lines in Fig. [5]). The P R (n) (not reported) display similar, 
though somehow smaller, deviations. Fig. [5^ and b refer 
to the middle of the system and r = 0, r = 0.5, respec- 
tively. Deviations decrease upon increasing the system 
sizes and are larger in the case of absorbing BC, as the 
particles have less time to "thermalize" . Fig. [5] c and d 
refer to the rightmost site, again for r = and r = 0.5. 
In this case there is no evidence of a convergence to- 
wards the equilibrium distribution, but this result is not 
earthshaking, as the density of particles decreases for in- 
creasing N. In other words, the deviations concern an 
asymptotically vanishing fraction of particles. 
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FIG. 5: (Color online) Distributions of left moving particles 
Pt{n) (broken lines) for i — N/2 and (a,b) and i = N (c,d). 
Dashed, dotted-dashed, and dotted curves correspond to TV = 
1000, 4000 and 16000, respectively; the solid line is P eq (n); 
see text for details. All curves are rescaled to span unit area. 

Conclusions - In this Letter we have determined the 
density profiles in the case of superdiffusive motion in 
an open, finite domain for arbitrary reflectivity r. Unex- 
pectedly, their shape depends on both a and r through 
the meniscus exponent /it, Eq. ([3]). This is a manifesta- 
tion of the long-range effects induced by superdiffusive 
behavior. The prediction of the Levy model accounts 
quantitatively for the shape of the temperature profile of 
an interacting oscillator chain displaying anomalous heat 
conductivity. It remains to clarify the relationship be- 
tween the type of BC in the thermal problem with that 
in the diffusive one, that appears to be far from trivial. 



Another open problem is the rigorous proof of the conjec- 
ture ([5]). We finally hope that the shape can be observed 
in some optical experiment. 

We thank J. Bertolotti for fruitful discussions. This 
work is part of the PRIN project 2008Y4W3CY. 
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